{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import sisl.viz\n",
    "import numpy as np\n",
    "from matplotlib import animation\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import HTML\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE**: Please complete exercise [TB 5](../TB_05/run.ipynb) before doing this exercise. \n",
    "\n",
    "In this example, you will learn how to find and calculate properties of *bound states* in NEGF calculations. In this example, the bound states are physically decoupled from the electrodes, hence *bound*. However, in more complicated calculations one might still encounter bound states, even if the atoms are neighbours. E.g. it may be that symmetries makes them *bound*.\n",
    "\n",
    "This example is very similar to [TB 5](../TB_05/run.ipynb), with the small change that some atoms are retained in the middle of the hole.\n",
    "You should again re-use the self-energies which, when built from scratch, require inverting, multiplying and adding matrices, roughly 10–20 times per $k$-point, per energy point, per electrode. \n",
    "For systems with large electrodes compared to the full device, this becomes more demanding than calculating the Green function for the system.  \n",
    "When there is periodicity in electrodes along the transverse semi-infinite direction (not along the transport direction) one can utilize Bloch's theorem to reduce the computational cost of calculating the self-energy.\n",
    "\n",
    "> In ***ANY*** calculation if you have periodicity, please ***USE*** it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(orthogonal=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the below lines are saving the electrode electronic structure *without* extending it 25 times (as was done in [TB_4](../TB_04/run.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_elec = sisl.Hamiltonian(graphene)\n",
    "H_elec.construct(([0.1, 1.43], [0., -2.7]))\n",
    "H_elec.write('ELEC.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = H_elec.tile(25, axis=0).tile(15, axis=1)\n",
    "\n",
    "# Get center of the atomic positions\n",
    "center = H.geometry.center(what=\"xyz\")\n",
    "# size of bound-state-flake\n",
    "bound_R = 5.\n",
    "\n",
    "_, ring_idx = H.geometry.close(\n",
    "    center,\n",
    "    R=[bound_R, 10.]\n",
    ")\n",
    "H = H.remove(ring_idx)\n",
    "\n",
    "dangling = [ia for ia in H.geometry.close(center, R=14.)\n",
    "                if len(H.edges(ia)) < 3]\n",
    "H = H.remove(dangling)\n",
    "edge = [ia for ia in H.geometry.close(center, R=14.)\n",
    "         if len(H.edges(ia)) < 4]\n",
    "edge = np.array(edge)\n",
    "\n",
    "# Pretty-print the list of atoms on the edges\n",
    "bound_idx = H.geometry.close(center, R=bound_R)\n",
    "# Make sure the edge atoms are only on the hole\n",
    "edge = np.setdiff1d(edge, bound_idx)\n",
    "print(\"All bound-state atoms: \", sisl.utils.list2str(bound_idx + 1))\n",
    "print(\"All edge atoms: \", sisl.utils.list2str(edge + 1))\n",
    "\n",
    "H.geometry.write('device.xyz')\n",
    "H.write('DEVICE.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "Instead of analysing the same thing as in [TB 5](../TB_05/run.ipynb) you should perform the following actions to explore the available data-analysis capabilities of TBtrans. *Remember*: always use Bloch's theorem when applicable!\n",
    "\n",
    "*HINT*: please copy as much as you like from example 04/05 to simplify the following tasks.\n",
    "\n",
    "1. Read in the resulting file into a variable called `tbt`.\n",
    "2. In the following, we will concentrate on *only* looking at $\\Gamma$-point related quantities. I.e. all quantities should only be plotted for this $k$-point.  \n",
    "   To extract information for one or more subset of points, you should look into the function\n",
    "       \n",
    "       help(tbt.kindex)\n",
    "   \n",
    "   which may be used to find a resulting $k$-point index in the result file.\n",
    "   \n",
    "3. Plot the transmission ($\\Gamma$-point only). To extract transmission for a subset $k$-point you should read the documentation of the specific function (*HINT*: `kavg` is the keyword you are looking for).\n",
    "   - Full transmission\n",
    "   - Bulk transmission. \n",
    "   How does it compare to example [TB 5](../TB_05/run.ipynb)? Should it be different?\n",
    "4. Plot the DOS with normalization according to the number of atoms ($\\Gamma$ only)  \n",
    "   Plot on the edge atoms, and compare with example [TB 5](../TB_05/run.ipynb)\n",
    "   - The Green function DOS (edge and bound atoms)\n",
    "   - The spectral DOS (edge and bound atoms). What is special about this quantity on the bound atoms?\n",
    "   - The bulk (electrode) DOS\n",
    "\n",
    "**TIME**: Do the calculation for various `TBT.Contour.Eta` values, in different folders. Then compare the same DOS as above. What changes? Why? Which DOS changes the most? The DOS on edge atoms or on bound atoms?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transmission"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Density-of-states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbt = sisl.get_sile('siesta.TBT.nc')\n",
    "# Easier manipulation of the geometry\n",
    "geom = tbt.geometry\n",
    "a_dev = tbt.a_dev # the indices where we have DOS\n",
    "# Extract the DOS, per orbital (hence sum=False)\n",
    "DOS = tbt.DOS(sum=False)\n",
    "# Normalize DOS for plotting (maximum size == 400)\n",
    "# This array has *all* energy points and orbitals\n",
    "DOS /= DOS.max() / 400\n",
    "a_xyz = geom.xyz[a_dev, :2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "fig = plt.figure(figsize=(12,4));\n",
    "ax = plt.axes();\n",
    "scatter = ax.scatter(a_xyz[:, 0], a_xyz[:, 1], 1);\n",
    "ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]');\n",
    "ax.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If this animation does not work, then don't spend time on it!\n",
    "def animate(i):\n",
    "    ax.set_title('Energy {:.3f} eV'.format(tbt.E[i]));\n",
    "    scatter.set_sizes(DOS[i]);\n",
    "    return scatter,\n",
    "anim = animation.FuncAnimation(fig, animate, frames=len(tbt.E), interval=100, repeat=False)\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learned lessons\n",
    "\n",
    "- Extraction of coupled elements via `.edges` for the Hamiltonian.\n",
    "- Manipulating the Hamiltonian *after* it has been created (*very* fast!)\n",
    "- Extract data only for single $k$-points. The lesson learned is also applicable for a subset of $k$-points.\n",
    "- Extraction of various physical quantities from the `*.TBT.nc` file\n",
    "- How to visualize bound states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
